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Abstract 

There is a whole range of emergent phenomena in non-equilibrium behaviors can be well de- 
scribed by a set of stochastic differential equations. Inspired by an insight gained during our study 
of robustness and stability in phage lambda genetic switch in modern biology, we found that there 
exists a classification of generic nonequilibrium processes: In the continuous description in terms 
of stochastic differential equations, there exists four dynamical elements: the potential function </>, 
the friction matrix S , the anti-symmetric matrix T, and the noise. The generic feature of absence 
of detailed balance is then precisely represented by T. For dynamical near a fixed point, whether or 
not it is stable or not, the stochastic dynamics is linear. A rather complete analysis has been car- 



ried out (Kwon, Ao, Thouless, cond-mat/0506280 ; PNAS, 102 (2005) 13029), referred to as SDS I. 



One important and persistent question is the existence of a potential function with nonlinear force 
and with multiplicative noise, with both nice local dynamical and global steady state properties. 
Here we demonstrate that a dynamical structure built into stochastic differential equation allows 
us to construct such a global optimization potential function. First, we provide the construction. 
One of most important ingredient is the generalized Einstein relation. We then present an ap- 
proximation scheme: The gradient expansion which turns every order into linear matrix equations. 
The consistent of such methodology with other known stochastic treatments will be discussed in 
next paper, SDS III; and the explicitly connection to statistical mechanics and thermodynamics 
will be discussed in a forthcoming paper, SDS IV. (The main results were published. Please cite 
the present paper as Potential in Stochastic Differential Equations: Novel Construction, P. Ao, J. 
Phys. A37 L25-L30 (2004). |http:/ / www.iop.org/EJ /abstract /0305-4470 /37/3 /L01/| ) 
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Let us consider an n component network whose dynamics is described by a set of stochastic 
differential equations 1 : 

Qtj = fjfat) +(j(lt,t). (1) 

The question is whether or not we can find a potential function from Eq.(l) which gives a 
global description of dynamics. Here qtj = dqtj/dt with j = 1, 2, n and the subscript t for 
q indicates that the state variable q is a function of time. The value of j th component is 
denoted by qj. The network state variable forms an n dimensional vector q T = (qi, q2, q n ) 
in the state space. Here the superscript r denotes the transpose. The state variable may 
be the values of particle coordinates in physics or the protein numbers in a signal trans- 
duction pathway, or any other possible quantities specifying the network. Let /j(q) be the 
deterministic nonlinear force on the j th component, which includes both the effects from 
other components and itself, and Q(q,t) the random force. For simplicity we will assume 
that fj is a smooth function explicitly independent of time. To be specific, the noise will be 
assumed to be Gaussian and white with the variance, 

(Ci(q t ,00(q*',O) = 2A;(q*W*- 0, (2) 

and zero mean, (Q) = 0. Here 8(t) is the Dirac delta function, and (...) indicates the 
average with respect to the dynamics of the stochastic force. By the physics and chemistry 
convention the semi-positive definite symmetric matrix D = {Dij} with i,j = 1,2, ...,n is 
the diffusion matrix. Eq.(2) also implies that, in situations where the temperature T can 
be defined, we have set k^T = 1 with k% the Boltzmann constant. We remark that if an 
average over the stochastic force (, a Wiener noise, is performed, Eq.(l) is reduced to the 
following equation in dynamical systems: 

(q*> = (f (%)> = f «*» • 

The last equality is due to the fact that at same time t, the noise and the state variables 
are independent of each other. It is equivalent to the fact that the noise can be switch 
out without affecting the deterministic force, a process demonstrated possible in physics in 
dealing with environmental effects 2 . A broad range of phenomena in both natural and social 
sciences has been described by such a deterministic equation 3 . 

Because of its importance and usefulness, repeated attempts have been made to construct 
a potential function^^^. The effort had, however, only limited success 8 . The usefulness of a 
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potential reemerges in the current study of dynamics of gene regulatory networks*^, which 
would again require its construction in complex network dynamics. It has been observed 
that the nonlinear dynamics is in general dissipative (tr(F) 7^ 0), asymmetric (Fji 7^ Fij), 
and stochastic (£ 7^ 0). Here the force matrix F is defined as 

Fij = dfi/dqj i,j = 1, ...,n . (3) 

and the trace is equal to the divergence of the force: tr(F) — d ■ f = J2]=idfj/dqj. The 
combination of those three features prevents any direct application of the insight from Hamil- 
tonian dynamics and has been the main obstacle preventing the potential construction. In 
fact, the asymmetry of dynamics has been characterized as the hallmark of the network 
in a state far from thermal equilibrium, and has been proclaimed that it makes the usual 
theoretical approach near thermal equilibrium unworkable 4 . It is the goal of this paper to 
report that we have, nevertheless, discovered a novel construction that can take care of those 
dynamical features and can give us a potential function. 

We state, the explicit construction will be given below, that there exists a unique decom- 
position such that Eq.(l) can be rewritten in the following form: 

[S(qt) + A(q t )]q t = -90(q t ) + £(q t , t) , (4) 

with the semi-positive definite symmetric matrix S'(qt), the anti-symmetric matrix A(qt), the 
single- valued scalar function 0(q*), and the stochastic force £(q$,t). Here d is the gradient 
operator in the state variable space. It is straightforward to verify that the semi-positive 
definite symmetric matrix term is 'dissipative': ifi.S{<\t)i\ t > 0; the anti-symmetric part does 
no 'work': q[y4(q t )qt = 0, therefore non-dissipative. Hence, it is natural to identify that the 
dissipation is represented by the semi-positive definite symmetric matrix S(q), the friction 
matrix, and the transverse force by the anti-symmetric matrix A(q), the transverse matrix. 
The scalar function 0(q) then acquires the meaning of potential energy. 

The decomposition from Eq.(l) to (4) may be called the ^-decomposition. However, 
without further constraint, Eq.(4) would be not unique. This may be illustrated by a simple 
counting: There are four apparent independent quantities in Eq.(4), while there are only 
two in Eq.(l). In order to have a unique form for Eq.(4), we may choose to impose the 
constraint on the stochastic force and the semi-positive definite symmetric matrix in the 
following manner: 

(a<lut)C(cit',t')} = 2S(q t )5(t-t') , (5) 
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and (£(qt, t)) = 0. We observe that this constraint is consistent with the Gaussian and white 
noise assumption for ( in Eq.(l). It may be called the stochasticity-dissipation relation. We 
further observe that that the forms of Eq.(4) and (5) strongly resemble those of dissipative 
dynamics in quantum mechanics when both dissipative and Berry phase exists^ 1 ^. The 
constrained 0- decomposition will be called the gauged 0-decomposition, which is indeed 
unique, as we will now demonstrate. 

We prove the existence and uniqueness of the gauged 0- decomposition from Eq.(l) to 
(4) by an explicit construction. Using Eq.(l) to eliminate the velocity q 4 in Eq.(4), we have 

[S(q t ) +A(q t )][f(q t ) +C(q*,*)] = -d<K<b) + £(q t , t) . 

Noticing that the dynamics of noise is independent of that of the state variables we require 
that both the deterministic force and the noise satisfying following two equations separately. 
For the deterministic force, this leads to 

[S(q)+A(q)]f(q) = -<90(q), (6) 

suggesting a 'rotation' from the force f to the gradient of the potential at each point in 
the state space. We have dropped the subscript t. For stochastic force, we have: 

[S(q)+A(q)]C(q,t)=£(q,t), (7) 

which shows the same 'rotation' between the stochastic forces. Here we have also dropped 
the subscript t for the state variable. Using Eq.(2) and (5), Eq.(7) implies 

[5(q) + A(q)]D(q)[5(q) - A(q)] = fif(q) , (8) 

which suggests a duality between Eq.(l) and (4): a large friction matrix implies a small 
diffusion matrix. It is a generalization of the Einstein relation- to finite transverse matrix 

Next we introduce an auxiliary matrix function 

G(q) = [5(q)+A(q)]- 1 . (9) 

Here the inversion ' — 1' is with respect to the matrix. Using the property of the potential 
function 0: d x <90 = 0, Eq.(6) leads to 

d x [G-H{q)\ = , (10) 
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which gives n(n — 1)/2 conditions to determine the nxn auxiliary matrix G. The generalized 
Einstein relation, Eq.(8), leads to the following equation 

G + G T = 2D, (11) 

which readily determines the symmetric part of the auxiliary matrix G, another n(n + l)/2 
conditions for G. Eq. (10) and (11) give the needed nxn conditions to completely determine 
G. Here we give a solution of G as a series in gradient expansion: 

oo 

G = D + J2^, (12) 

3=0 

with AGj = Y£i(-l) l [(F T ) l D3F~ l + {F T )~ l DjF% D = DF — F T D, = (D + 

AGj-i) {[d x (D- 1 + AGjl x )]f} {D - AG J _ 1 ). The zeroth order approximation to Eq.(10) 
is GF T — FG T = 0. A formal solution to this approximated equation together with Eq.(ll) 
has been constructed under a rather restrictive condition 1 -, and the explicit solution under 
a generic condition has been obtained in Ref. [14]. The gauged (^-decomposition is therefore 
uniquely determined: 

'0(q) = -Jcdq'-^qW)] 

< 5(q) = [G-\ q ) + (GT)-\ q )]/2 . (13) 
k A(q) = [G-\ q ) - (Cr)-\<0]/2 
The end and initial points of the integration contour C are q and qo respectively. During the 
construction a sufficient condition det(F) det(S , + y4) ^ is assumed, with exception at some 
isolated points. We remark on the special role played by the force matrix: If FD = DF T , 
i.e., D = 0. If in addition Dj>i = in this case, which can be realized if D is a constant, 
AGj = for all j. This means that G = D and the transverse matrix A = 0. Such a 
condition has been noticed in the linear case where both F and D are constant matrices, 
and named the integrability condition^. 

In many experimental studies of a complex network, a question is often asked on the 
distribution of the state variable after a transient period instead of focusing on the individual 
trajectory of the network. This implies that either there is an ensemble of identical networks 
or repetitive experiments are been carried out. From statistical mechanics, if viewing the 
potential function as an energy, a steady distribution function can be expected from Eq.(4): 

Po(q) = ^ exp{-0(q)} , (14) 
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with the partition function Z = J cf n q exp{— </>(q)}. This is a Boltzmann-Gibbs distribution 
for the state variable, and give the strongest manifestation of the usefulness of the potential 
function <fi. We remark that it is however not obvious that the steady state distribution, if 
exists, should be given by Eq.(14). In the following we give a heuristic demonstration that 
Eq.(14) is indeed a right steady distribution for the network as the steady state solution of 
the corresponding Fokker-Planck equation. 

The connection between the standard stochastic differential equation, Eq.(l), and Fokker- 
Planck equation is necessarily ambiguous for the generic nonlinear exemplified by 
the Ito-Stratonovich dilemma^^. We attribute this lack of definiteness to the asymptotic 
nature of the connection in which a procedure must be explicitly defined: Different proce- 
dures will in general lead to different results. Here we present another procedure which may 
be natural from the theoretical physics point of view. Our starting point will be Eq.(4), not 
Eq.(l) from which most previous derivations started. 

The existence of both the deterministic force and the stochastic force in Eq. (4) suggests 
that there are two well separated time scales in the network: microscopic time scale for 
the description of the stochastic force and macroscopic time scale for the network motion. 
The former time scale is much smaller than the latter. This separation of time scales 
further suggests that the macroscopic motion of the network has an inertial: It cannot 
response instantaneously to the microscopic motion. To capture this feature, we introduce a 
small inertial "mass" m and a kinetic momentum vector p for the network. The dynamical 
equation for the network now takes the form: 



is the extension of Eq.(4). We note that there is no dependent of friction matrix and the 
stochastic force on the kinetic momentum, therefore no Ito-Stratonovich dilemma in the 
connection between the stochastic differential equation and the dynamical equation for the 
distribution function. The Fokker-Planck equation in this enlarged state space, the Klein- 
Kramers equation, can be immediately obtained 1 : 



q t = pt/m 



(15) 



defines the kinetic momentum, and 




(16) 




(17) 
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Here f = pA/m — <9 q </>, and t, q, and p are independent variables. The steady distribution 
can be found as^ 

p(q, p) = exp{-[p 2 /2m + 0(q)]}/Z , (18) 

with Z = J d n qd n p exp{ — [p 2 /2m + 0(q)]} the partition function. There is an explicit 
separation of state variable and its kinetic momentum in Eq.(18). The zero "mass" limit 
can then be taken with no effect on the state variable distribution. This confirms that 
Eq.(14), the Boltzmann-Gibbs distribution expected from Eq.(4), is the right choice under 
this procedure. 

To conclude, we point out a major difference between the present construction of the 
potential and those in literature such as represented by the Graham-Haken construction^: 
The present construction is based on a structure built into stochastic differential equation. 
There is no explicit use of Fokker-Planck equation. Therefore there is no need to make 
assumption on the distribution function in the limit time goes to infinite as assumed in the 
Graham-Haken construction. In particular, the potential in the present paper can be time 
dependent. 

Critical discussions with N. Goldenfeld, A.J. Leggett, C. Kwon, J. Ross, D.J. Thouless, 
L. Yin, and X.-M. Zhu are highly appreciated. This work was supported in part by a USA 
NIH grant under the number HG002894-01. 
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